function SS = skill_score(modeldata, obsdata, do_SS_breakdown)
%
% this computes the skill score based on two time series
%
% SS = skill_score(modeldata, obsdata, do_SS_breakdown)
%  as SS = 1 - (1/var(o))*(1/N)*sum([m - o]^2) 
%  
% do_SS = 1 to get terms in SS (like CC2, etc.)
%
% where var is varianece, N is length of time series
% and m = modeled, o = observed,
% (Murphy 1988)
% from DAS, updated by SNG May 2011
%************************************************

if(nargin<3)
    do_SS_breakdown = 0; %only report SS, not terms in SS
end


m = modeldata(:);
o = obsdata(:);
good = ~isnan(m) & ~isnan(o);
m = m(good);
o = o(good);

if(length(m)~=length(o))
    error('Vectors must be the same length!');
end

varO = nanvar(o);  %obs variance
varM = nanvar(m);  %obs variance
N = length(o);

MSE_top = nanmean((m - o).^2);
MSE_bot = nanmean((o - nanmean(o)).^2);

if(~do_SS_breakdown)
   SS = 1 - MSE_top ./ MSE_bot;
elseif(do_SS_breakdown)
   SS.SS = 1 - MSE_top ./ MSE_bot;
   %cc = corr_coef(m,o);
   r = corrcoef(m,o); cc = r(1,2);
   meanO = nanmean(o); meanM = nanmean(m);
   % terms in SS, CC^2, variance term, mean term
   SS.CC2 = cc.^2;
   SS.var = (cc - sqrt(varM)./sqrt(varO)).^2; 
   SS.mean = ((meanM-meanO)./sqrt(varO)).^2;
   % where SS = CC2 - var - mean 
end

% gives skill score on whether model is better than mean of observations
% (0 is no better, >0.2 is good, <0 is worse than obs.)





